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Many mesh refinement simulations currently performed in numerical relativity counteract insta- 
bilities near the outer boundary of the simulation domain either by changes to the mesh refinement 
scheme or by changes to the gauge condition. We point out that the BSSN Gamma Driver gauge 
condition introduces a time step size limitation in a similar manner as a CFL condition, but which 
is independent of the spatial resolution. We give a didactic explanation of this issue, show why 
especially mesh refinement simulations suffer from it, and point to a simple remedy. 



I. INTRODUCTION 
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We assume that the audience is familiar with the con- 
cept of a Courant-Friedrichs-Lewy (CFL) condition 
Loosely speaking, the CFL condition states: When a par- 
tial differential equation, for example the wave equation 



Am 



(1) 



is integrated numerically, then the time step size St is 
limited by the spatial resolution Sx and the maximum 
propagation speed c by 



St < Q 



Sx 



(2) 



Here Q is a constant of order 1 that depends on the time 
integration method (and details of the spatial discreti- 
sation). Choosing a time step size larger than this is 
unstable and must therefore be avoided. (There are time 
integration methods that do not have such a stability 
limit, but these are expensive and not commonly used in 
numerical relativity, so we will ignore them here.) 



II. EXAMPLE: EXPONENTIAL DECAY 

In real- world equations, there are also other restric- 
tions which limit the time step size, and which may be 
independent of the spatial resolution. One simple exam- 
ple for this is the exponential decay 



d t u 



-A u 



(3) 



where A > is the decay constant. Note that this equa- 
tion is an ordinary differential equation, as there are no 
spatial derivatives. The solutions of ^ are given by 



ult) = A cxp{-Ai} 



(4) 



with amplitude A. 
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The decay constant A has dimension 1/T. The time 
step size is limited by 



St < Q' 



(5) 



where Q' is a constant of order 1 that depends on the 
time integration method. Choosing a time step size larger 
than this is unstable and must therefore be avoided. (As 
with the CFL criterion, there are time integration meth- 
ods that do not have such a stability limit.) 

As an example, let us consider the forward Eulcr 
scheme with a step size St. This leads to the discrete 
time evolution equation 
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or 



= -Ai 
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This system is unstable e.g. if |u ri+1 | > 
also other definitions of stability), or if 

\l-StX\ > 1 , 



(6) 
(7) 

(there are 
(8) 



which is the case for St > 2/ A (and also for St < 0). 
In this case, the solution oscillates between positive and 
negative values with an exponentially growing amplitude. 



III. GAMMA DRIVER 

The BSSN Q Gamma Driver condition is a time evo- 
lution equation for the shift vector /3 1 , given by (see e.g. 
(43) in 1) 



(9) 



There exist variations of the Gamma Driver condition, 
but the fundamental form of the equation remains the 
same. The term FdtT 1 contains second spatial deriva- 
tives of the shift /3 l and renders this a hyperbolic, wave- 
type equation for the shift. The parameter rj > is 
a damping parameter, very similar to A in {3| above. 
It drives 9 t /3 l to zero, so that the shift (3 1 will tend to 
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a constant in stationary spacetimcs. (This makes this 
a symmetry-seeking gauge condition, since d t will then 
tend to the corresponding Killing vector.) 

Let us now consider a simple spacetime which is spa- 
tially homogeneous, i.e. where all spatial derivatives van- 
ish. In this case (see e.g. (40) in |3|), dtT 1 = 0, and only 
the damped oscillator equation 



(10) 



remains. As we have seen above, solving this equation nu- 
merically still imposes a time step size limit, even though 
there is no length scale introduced by the spatial dis- 
cretisation, so the spatial resolution can be chosen to be 
arbitrarily large; there is therefore no CFL limit. This 
demonstrates that the damping time scale set by the pa- 
rameter rj introduces a resolution-independent time step 
size limit. 

This instability was e.g. reported in below (13) 
there, without explaining its cause. The authors state 
that the choice 77 = 2 is unstable near the outer bound- 
ary, and they therefore choose rj = 1 instead. Decreasing 
77 by a factor of 2 increases the time step size limit cor- 
respondingly. 

The explanation presented above was first brought 
forth by Carsten Gundlach [f| and Ian Hawke @. To 
our knowledge, it has not yet been discussed in the liter- 
ature elsewhere. 

Harmonic formulations of the Einstein equations have 
driver parameters similar to the BSSN Gamma Driver pa- 
rameter T). Spatially varying parameters were introduced 
in harmonic formulations to simplify the gauge dynamics 
in the wave extraction zone far away from the origin (see 
e.g. (8) in @). Q uses a harmonic formulation with mesh 
refinement, and describes using this spatial dependence 
also to avoid time stepping instabilities (see (45) there). 



IV. MESH REFINEMENT 

When using mesh refinement to study compact objects, 
such as black holes, neutron stars, or binary systems of 
these, one generally uses a grid structure that has a fine 
resolution near the centre and successively coarser reso- 
lutions further away from the centre. With full Berger- 
Oliger AMR that uses sub-cycling in time, the CFL fac- 
tors on all refinement levels are the same, and thus the 
time step sizes increase as one moves away from the cen- 
tre. This makes it possible that the time step size on the 
coarsest grids does not satisfy the stability condition for 
the Gamma Driver damping parameter 77 any more. 

One solution to this problem is to omit sub-cycling in 
time for the coarsest grids by choosing the same time 
step size for some of the coarsest grids. This was first 
advocated by 0, although it was introduced there to 
allow large shift vectors near the outer boundary as nec- 
essary for a co-rotating coordinate system. It was later 
used in [l(| (see section IV there) to avoid an instabil- 
ity near the outer boundary, although the instability is 



there not attributed to the Gamma Driver. Omitting 
sub-cycling in time on the coarsest grids often increases 
the computational cost only marginally, since most of the 
computation time is spent on the finest levels. 

Another solution is to choose a spatially varying pa- 
rameter 77, e.g. based on the coordinate radius and mim- 
icking the temporal resolution of the grid structure, 
which may grow linearly with the radius. This follows 
the interpretation of rj setting the damping timescale, 
which must not be larger than the timescale set by the 
time discretisation. 

One possible spatially varying definition for 77 could be 
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where r is the coordinate distance from the centre of the 
black hole. The parameter R defines a transition radius 
between an inner region, where 77 is approximately equal 
to 77*, and an outer region, where 77 gradually decreases 
to zero. This definition is simple, smooth, and diffcr- 
cntiablc, and mimics a "typical" mesh refinement setup, 
where the resolution h grows approximately linearly with 
the radius r. 

Another, simpler definition for 77 (which is not smooth 
- but smoothness is not necessary; 77 could even be dis- 
continuous) is 



77(7-) 



1 for r < R (near the origin) , . 
-p for r > R (far away) ' 



which is e.g. implemented in the McLachlan code [111 ]. 

If there are multiple black holes, possibly with differing 
resolution requirements, then prescriptions such as (|11[) 
or (fT2"j) need to be suitably generalised, e.g. via 
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77(7-) 771(7-1) 772(7-2) 



(13) 



where 771 and 772 are the contributions from the individ- 
ual black holes, with n and 7-2 the distances to their 
centres. This form of (fT3l is motivated by the dimension 
of 77, which is 1/M, so that two superposed black holes 
of masses mi and 7712 lead to the same definition of 77 as 
a single black hole with mass mi + 777.2. 

Another prescription for a spatially varying 77 has been 
suggested in [l2|. In this prescription, 77 depends on the 
determinant of the three-metric, and it thus takes the 
masses of the black hole(s) automatically into account. 
This prescription is motivated by binary systems of black 
holes with unequal masses, where 77 near the individual 
black holes should be adapted to the individual black 
holes' masses, and it may be more suitable to use this 
instead of JT3|). 

There can be other limitations of the time step size 
near the outer boundary, coming e.g. from the boundary 
condition itself. In particular, radiative boundary condi- 
tions impose a CFL limit that may be stricter than the 
CFL condition from the time evolution equations in the 
interior. 
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